Computational investigation of IP3 diffusion

Inositol 1,4,5-trisphosphate (IP3) plays a key role in calcium signaling. After stimulation, it diffuses from the plasma membrane where it is produced to the endoplasmic reticulum where its receptors are localized. Based on in vitro measurements, IP3 was long thought to be a global messenger characterized by a diffusion coefficient of ~ 280 μm2s−1. However, in vivo observations revealed that this value does not match with the timing of localized Ca2+ increases induced by the confined release of a non-metabolizable IP3 analog. A theoretical analysis of these data concluded that in intact cells diffusion of IP3 is strongly hindered, leading to a 30-fold reduction of the diffusion coefficient. Here, we performed a new computational analysis of the same observations using a stochastic model of Ca2+ puffs. Our simulations concluded that the value of the effective IP3 diffusion coefficient is close to 100 μm2s−1. Such moderate reduction with respect to in vitro estimations quantitatively agrees with a buffering effect by non-fully bound inactive IP3 receptors. The model also reveals that IP3 spreading is not much affected by the endoplasmic reticulum, which represents an obstacle to the free displacement of molecules, but can be significantly increased in cells displaying elongated, 1-dimensional like geometries.

where D represents the IP 3 diffusion coefficient in the absence of buffers, S T the concentration of IP 3 R monomers, and K D the equilibrium dissociation constant of IP 3 from its receptor. S T is cell dependent, in the range of 80 nM to 2 μM 13, [16][17][18] with an estimated value of 542 nM in SH-SY5Y cells 18 . Thus, the buffering effect described by Eq. 1 is expected to induce a reduction of the effective diffusion coefficient of ~ 2.25 at [IP 3 ] = K D . In Dickinson et al. 11 experiments, [IP 3 ] was in the Ca 2+ oscillatory range. Because IP 3 concentrations are typically of the order or lower than the K D of the IP 3 R in the oscillatory range 19 , one can estimate that the effective diffusion coefficient of IP 3 considering IP 3 buffering by not fully bound receptors (D I ) should be ~ 100 μm 2 s -1 (Fig. 1). Values around 5-10 μm 2 s -1 would correspond to situations in which IP 3 concentrations are lower than K D /2 and IP 3 R concentrations at least higher than 1 μM, in which IP 3 binding to its receptors exceeds what occurs in SH-SY5Y cells. Such discrepancy led Dickinson et al. 11 to hypothesize the existence of two distinct populations of IP 3 Rs, with the most active exhibiting different binding and/or gating properties. These different populations were proposed to be the molecular bases of the observed two modes of Ca 2+ release: one early pulsatile Ca 2+ activity at puff sites and one later, spatially more diffuse Ca 2+ liberation 20 .
Besides, a value of D I lower than 10 μm 2 s −1 raises questions about some observations related to Ca 2+ waves. For example, in ascidian eggs that have a radius between 25 and 38 μm, Ca 2+ waves propagate across the entire fertilized egg in less than 10 s while IP 3 is only synthesized at the plasma membrane 21 . The question is even more critical for intercellular Ca 2+ waves, which in many cases rely on the propagation of IP 3 through gap junctions 22 . Intracellular diffusion of IP 3 with a coefficient ≤ 10 μm 2 s −1 could not account for waves propagating at rates ≥ 10 μms −1 among a cell population 23,24 . More generally, an accurate quantitative description of IP 3 diffusion is necessary for the numerous modelling studies devoted to the physiological responses induced by IP 3 -controlled Ca 2+ signaling, such as, for example, the Ca 2+ -and IP 3 -regulated nitric oxide production in neurons 25 , the metabolism of β-amyloids during the development of Alzheimer's disease 26 or salivary secretion by acinar cells 27 . IP 3 spreading from its location of synthesis to its receptors also plays a key role in cardiac cells 28,29 .
The aim of the present study is to re-investigate the characteristics of IP 3 diffusion using a stochastic model that explicitly simulates Ca 2+ puff dynamics and allows a realistic computational description of the experiments performed by Dickinson et al. 11 . Based on the observations of puff latencies performed by these authors, we propose a new computational treatment to infer the value of the effective diffusion coefficient of IP 3 , leading to different conclusions. The reasons for the different outcomes between the two studies are analyzed in the discussion. We also investigate how cell geometry and the presence of the ER membranes, which provide an obstacle to the free displacement of molecules, affect the rates at which IP 3 diffuses in a cell-like environment.

Results
Methodology. We performed explicit stochastic simulations of Ca 2+14 releasing activities of puff sites located at different distances from an IP 3 source, in order to directly simulate the experimental protocol that was used to estimate the IP 3 diffusion coefficient in vivo 11 . The mathematical model is based on a previously proposed fully stochastic description of the Ca 2+ exchanges between the ER and the cytosol via IP 3 Rs, SERCA pumps and a leak from the ER 9 . In this work, SERCA pumps, Ca 2+ leakage and Ca 2+ diffusion are described deterministically. The effective value of the Ca 2+ diffusion coefficient taking Ca 2+ buffering into account 10 is considered, i.e. 40 μm 2 s −1 . To simulate Ca 2+ puffs, the release of Ca 2+ via the IP 3 Rs is described stochastically, using the Gillespie's algorithm 30 . Each cluster of IP 3 Rs (Supplemental Figure S1) is described as a whole and can be in four states: one open (O), one closed (C) and two different inhibited ones (I 1 and I 2 ). These states describe the global behavior of a cluster composed of close-by IP 3 Rs. The phenomenological model was shown to reproduce experimentally observed statistical properties of Ca 2+ puffs 31 and to describe the passage from localized puffs to global Ca 2+ spikes when clusters are effectively coupled by Ca 2+ diffusion 9 .
To evaluate the effective diffusion coefficient of IP 3 , this model is extended to take the dependence of puff occurrence on IP 3 concentration into account. Thus, in the Gillespie's simulations, the propensity of transition of the cluster from the closed (C) to the open (O) state now writes: where k CO stands for the rate constant that characterizes the passage of the cluster from the closed to the open state, N Ca is the number of cytosolic Ca 2+ ions and Ω is the extensivity parameter. N IP and K IP represent the number of IP 3 molecules and the IP 3 dissociation constant of the IP 3 R (multiplied by Ω), respectively. Although the model does not simulate individual IP 3 receptors, we assumed that puff firing probability is related to the probability of one tetrameric IP 3 R to be fully bound to IP 3 13 . As shown below, Eq. 2 indeed allows to reproduce the exponential dependence of mean first puff latency on IP 3 concentration reported experimentally 11,32 .
The evolution of IP 3 concentration is described deterministically and two different protocols of photorelease of caged IP 3 are simulated. The first one corresponds to a spot photorelease of caged IP 3 in a small region of the cell: In the second one, called "distributed photorelease", IP 3 is liberated at different spots to get an increase in [IP 3 ] that is nearly spatially homogeneous in the whole cell. In this case: In Eqs. (3)(4)(5)(6), D I stands for the effective IP 3 diffusion coefficient, θ for the rate of IP 3 release upon laser flash and H for Heaviside function. This function equals 1 if t ≤ t ph and 0 otherwise. Considering that θ is determined by the intensity of the flash in the spot photorelease case (Eq. 3), a ten times smaller value is used to simulate distributed photorelease (Eq. 5) because the same total amount of IP 3 is distributed among the ten IP 3 releasing sites. Basal IP 3 concentration is set to 50 nM 33 .
A detailed description of the algorithm is provided in the Supplementary Information (Sect. 1), together with a Table listing the propensities used for the stochastic part of the algorithm and the evolution equations used in the deterministic part (Table Supplement 1). Code are available at https:// github. com/ Rober toOrn elasG uevara/ ca2-puffs/ tree/ main. Except for the process describing IP 3 dynamics listed here above, the values of parameters are taken from Voorsluijs et al. 9 and are listed in Table Supplement 2. Validation of the model. Low concentrations of IP 3 typically evoke local Ca 2+ signals known as Ca 2+ puffs 12,34 . As the concentration of IP 3 increases, Ca 2+ puffs become more frequent and transform into Ca 2+ waves spreading regeneratively across the cell. These repetitive waves, also known as Ca 2+ spikes, are more regular than puffs and are often referred to as Ca 2+ oscillations. Their stochastic origin is visible by the linear relation between the variance on the interspike interval and the mean interspike interval itself 8,35 . In experiments, Ca 2+ waves propagation can be hindered by loading the cells with the slow Ca 2+ buffer EGTA, which reduces communication between puff sites 36 . At high EGTA concentrations, each cluster practically behaves as an independent entity and thus generates Ca 2+ puffs on a large range of [IP 3 ].
Such behavior is well reproduced by the stochastic model, as seen by simulations in a simplified 2D geometry ( Fig. 2A-D). Four different conditions were used: low and intermediate IP 3 concentration, with and without coupling between clusters by Ca 2+ diffusion. In all cases, [IP 3 ] is constant in time and space. To model puff dynamics when clusters are uncoupled (corresponding to the presence of EGTA), a single cluster was simulated and the local Ca 2+ concentration averaged in 1fL around the cluster was monitored. The resulting dynamics of [Ca 2+ ] ( Fig. 2A,B) corresponds to Ca 2+ puffs, with statistical properties in good agreement with observations and an  . For the 2 panels, one cluster, located at the centre of a 5 × 5 μm 2 2D system is considered. For panels (C and D), a larger 5 × 10 μm 2 system containing 10 clusters is simulated. All values of parameters are the same as in A and B and listed in Table S2. Time series show the evolution of Ca 2+ concentration averaged on the whole system. Panels (E and F) show the spatio-temporal evolutions of Ca 2+ puffs and spikes, respectively. In the two cases, the system is 50 × 10 μm 2 large and 10% of the total surface is occupied by clusters (200 clusters). www.nature.com/scientificreports/ an isolated cluster of IP 3 Rs. To model spikes dynamics when clusters are coupled by diffusion, 10 clusters of IP 3 Rs were randomly distributed in a 5 μm × 5 μm system and the global [Ca 2+ ] averaged over the whole system was monitored (Fig. 2C,D). The mean interspike interval and the coefficient of variation (CV) decrease when increasing [IP 3 ], as reported previously 35  To simulate the spatio-temporal profiles of puffs and spikes, simulations were performed in a 50 μm × 10 μm system, while keeping the same average density of clusters (Fig. 2E,F). While a low, homogeneous concentration of IP 3 induces random puff activity, a large bolus of IP 3 released at one extremity of the cell initiates a global Ca 2+ spike that propagates as a Ca 2+ wave. Upon a spatially uniform increase in IP 3 concentration, latencies of first puff occurrence are exponentially distributed 12,32 . Such distributions, obtained by model simulations, are shown in Fig. 3A-C, together with their means and associated standard error of the mean (SEM) values. First puff latencies are exponentially distributed with an average that decreases with the IP 3 concentration, consistent with observations in SH-SY5Y cells 11 . The characteristic decrease times of the exponential distributions (τ) also decrease with [IP 3 ]. These results were used to infer the values of the rate of release of caged IP 3 in response to the flash, i.e. the value of parameter θ in Eqs. 3 and 5. Starting from the relation between mean first puff latency and flash duration reported by Dickinson et al. 11 for the distributed photorelease of IP 3 -and replotted in Fig. 3D, we seek the values of IP 3 concentrations that, in www.nature.com/scientificreports/ the simulations, gave the same mean latencies as those reported experimentally using Fig. 3C. Next, we numerically evaluated the value of θ allowing to reach these spatially uniform concentrations of IP 3 when assuming 10 photorelease spots of durations equal to 0.1, 0.2 and 0.5 s respectively, as in the experiments. Good agreement between mean latencies and flash durations was found for θ = 600 μM/s (Fig. 3D). When caged IP 3 is photoreleased at one extremity of the cell, puffs begin on average after longer latencies at greater distance from the spot 11 . This reflects the time taken for IP 3 to diffuse up to the cluster, and IP 3 dilution. Dickinson et al. 11 used mean first puff latencies at different distances from the photorelease spot to estimate the effective diffusion coefficient of IP 3 . In addition to mean first puff latencies, the time interval between the flash and the observation of the first puff at a given distance from the release site is also a relevant quantity. We called this latency "minimal first puff latency". It can be regarded as the time lapse during which all clusters located at a given distance from the flash site remain silent, indicating that IP 3 diffusion up to this point has been negligible. As visible in the distribution of puff latencies upon a global IP 3 increase, the minimal first puff latency is shorter than 0.5 s as soon as [IP 3 ] exceeds 80 nM (Fig. 3C, blue dashed line). In agreement with this, the minimal first puff latency is not affected by the duration of the flash in the case of distributed photorelease of IP 3 11 . If the source of IP 3 is spatially restricted, minimal first puff latency provides a reliable indication of the time at which IP 3 has started to increase at a given location since IP 3 diffusion can be viewed as a deterministic process and since the number of puff sites analyzed is large enough. Because puff activity is a stochastic process that also depends on the Ca 2+ concentration around the cluster, minimal first puff latencies in principle overestimate the time required for IP 3 to spread up to a given location. However, given that minimal first puffs are shorter than 0.5 s at any IP 3 concentration, this delay is anyway small compared to the time needed for IP 3 to spread over distances larger than several microns. Moreover, as discussed below, the minimal first puff latency is not significantly influenced by the possible incomplete Ca 2+ buffering by EGTA, which could accelerate puff triggering. Indeed, first puffs arise in conditions where most, if not all, nearby clusters are inactive.
Effective diffusion coefficient of IP 3 . We investigated the relation between first puff latency and distance from spot photorelease by performing independent simulations for twelve clusters, each one located at a different distance from the photorelease spot in an ellipse-shaped 2D cell. Such a configuration allows to reproduce the absence of communication between the cluster sites that results from the addition of EGTA. Binding of IP 3 to IP 3 Rs is taken into account in the value of the effective IP 3 diffusion coefficient, D I (see Eq. 1). The ellipse shape was chosen instead of the rectangle used in the previous figures to avoid artefactual effects in the corners. Each simulation was run up to the opening of the cluster and this time was then monitored. The same simulation was performed 50 times for each distance between the spot and the cluster. The averages of the latencies of the first puffs correspond to the mean latencies. For the minimal first puff latencies, we divided the 50 simulations in five groups of ten. In each group, the minimal value was spotted. The average of these minimal values is defined as the minimal first puff latency. As shown in Fig. 4 (upper row), minimal first puffs latencies (black stars) simulated with an effective IP 3 diffusion coefficient of 100 μm 2 s −1 -which corresponds to the expected value given the IP 3 buffering capacity of the cytosol of SH-SY5Y cells, as discussed in the Introduction-are in agreement with experimental observations of Dickinson et al. 11 . Thus, the time taken by IP 3 to diffuse from the photorelease spot to the locations of the clusters is such that the relation between the distance from the flash spot and the duration of the period of inactivity are well reproduced. Agreement between simulated and observed mean first puff latencies is limited to the case of the 0.1 s flash, i.e. the lowest IP 3 concentration. For larger IP 3 concentrations (0.2 and 0.5 s flash), we reasoned that the discrepancy between the idealized simulations of cells containing one single cluster and experiments may be due to a slight stimulation of puff activity by Ca 2+ in the experiments, which is not totally buffered by EGTA. At 5 μM EGTA, the distance on which free Ca 2+ ions diffuse before being captured by EGTA is indeed between 2.1 and 3.8 μm (see Sect. 2 of Supplementary Information). As these distances are of the order of the mean distance between clusters in the SH-SY5Y cells, there is some Ca 2+ stimulation of puff activity in the distributed photorelease protocol. Because we only model one cluster at a time, this effect does not occur in our simulations. Thus, the value θ = 600 μMs −1 deduced in the previous section overestimates the IP 3 release rate and must be seen as an upper limit.
In line with this assumption, a lower value of parameter θ (250 μMs −1 ) representing the rate of IP 3 release due to the flash (Eqs. 3 and 5) allows to get a good agreement both for the mean and the minimal first puff latencies (Fig. 4, medium row). This value of θ also allows us to reproduce the experimentally observed dependency of the mean latency on the flash duration when considering, in the simulations, that the Ca 2+ level is increased at cluster location during puff activity upon distributed IP 3 photorelease (Supplemental Figure S2). It should be stressed that the 250 μMs −1 value for θ was found to provide a good fit with observations but remains arbitrary since it cannot be computed on the basis of a detailed quantitative knowledge of the level of communication between the clusters via Ca 2+ . However, the third row of Fig. 4 shows that agreement with observations cannot be obtained when keeping the large value of θ (600 μMs −1 ) and decreasing the value of D I . In this case, despite the agreement for mean latencies, there is a systematic and important overestimation of the minimal first puff latencies. Such disagreement cannot be ascribed to a reduced number of experimental observations because a larger number of data is only expected to further decrease the shortest possible latency at a given distance. In contrast, it reveals that the low value of D I is unable to simulate the quite rapid, moderate IP 3 increase away from the spot. The larger concentration of IP 3 due to the overestimated θ explains why mean latencies agree with observations, as they compensate in average for the slow diffusion. Simulations also predict that with θ = 250 μMs −1 best agreement is obtained with D I = 100 μm 2 s −1 , as compared to smaller or larger values (Fig. 5). Altogether, simulations are in good agreement with experimental observations concerning mean and minimal first puff latencies with the value of the IP 3 diffusion coefficient taking IP 3  www.nature.com/scientificreports/ latter value also allowed to reproduce the relation between puff latencies and distance from the flash measured in COS-7 cells (Supplemental Figure S3).

Influence of cell shape.
Up to this point, simulations have been performed in 2D, rectangular or ellipseshaped systems. However, the 3D character and the specific geometry of the cells are expected to influence the rate at which IP 3 propagates into the cytoplasm through diffusion. To address this question using computational simulations (Fig. 6), we first considered an ellipsoid (50 × 10 × 5 µm) in which 20 cluster sites are located randomly at a distance shorter than 0.1 μm from the plasma membrane 37 . IP 3 is assumed to be released in a 0.25 μm radius sphere located at the left extremity of the cell. The rate of IP 3 increase (parameter θ in Eqs. 3 and 5) was adapted according to the change in the cytoplasmic volume. Two seconds after the simulated flash, a gradient of IP 3 is established (Fig. 6). Consequently, the dynamics of puff activity is rather different depending on the distance from the flash (compare the Ca 2+ time series in Fig. 6). Agreement between minimal and mean first puff latencies computed in simulations with an effective diffusion coefficient of IP 3 equal to 100 μm 2 s −1 and observations of Dickinson et al. 11 is slightly improved in this 3D configuration as compared to the 2D situation (Supplemental Figure S4). The shape of the cell is expected to influence the rate of IP 3 spreading. We investigated this possible influence by simulating IP 3 diffusion in a two-dimensional representation of an astrocyte (Fig. 7). In response to a release of   Fig. 4. The rate of localized IP 3 photorelease, θ, is taken equal to 2500 μMs −1 , which allows to obtain the same spatiotemporal profile of IP 3 increase as the 250 μMs −1 value for the 2D case, because of changes in volume. www.nature.com/scientificreports/ IP 3 at the intersection between the cell body and one process, Ca 2+ puffs occur sooner in the cell process than in the body. This is related to a larger rate of IP 3 spreading in elongated structures where dilution is much reduced. For example, 2 s after the simulated flash, more elevated IP 3 concentrations are seen in the process than in the cell body (Fig. 7D). Thus, at the same time and the same distance from the site of IP 3 release, IP 3 reaches higher concentrations in the process and is able to propagate on longer distances (Supplemental Video S5).
In conclusion, the value inferred for the effective IP 3 diffusion coefficient is not significantly affected by considering a 3D system instead of a 2D one. On the other hand, cell shape can considerably affect spreading, with elongated geometries increasing concentration gradients and thereby favoring fast diffusion.
Realistic ER geometry. The ER consists in a network of tubules and flattened sacs, resulting in a complex shape. Since the IP 3 receptors are located within its membrane, spreading of IP 3 to the receptors may be affected by the ER structure. Taking advantage of the great flexibility of Comsol Multiphysics, we performed simulations in which the ER structure was explicitly considered. Thus, we simulated the geometry of a SH-SY5Y cell, in which we inserted an ER whose shape was largely inspired from 2D images of ER in DC-3F cells 39 . We first compared the puff latencies simulated in the absence and in the presence of ER. The values of the rates of IP 3 release were kept the same as in the simulations above (Fig. 4), using a scaling factor to take into account the changes in the cytosolic volume. Puff latencies were in average not much affected by the presence of the ER (Fig. 8). However, the time evolution of the corresponding IP 3 profiles indicates that spreading is locally affected by the extent and the shape of the accessible portions of cytosol between the IP 3 release point and a given puff location. This is most easily visualized in an idealized, ellipse-shaped geometry. In this case, the IP 3 profiles along a fictive line located at half the cell length and perpendicular to the gradient are indeed slightly different when IP 3 is released www.nature.com/scientificreports/ at the right or the left extremity of the cell (Fig. 9A-C). Thus, IP 3 R clusters locally experience somewhat different IP 3 increases depending on the ER shape, which represents an additional source of randomness in puff activity. When compared to simulations where the ER structure is not considered, average IP 3 increases are in average either a bit slower or similar (Fig. 9D). However, locally, the presence of the ER can also allow for a faster IP 3 increase (Fig. 9B).
In conclusion, simulations indicate that the ER does not much affect average spreading of IP 3 in the cytosol. Although it constitutes a physical barrier that can locally slow down spreading, it also creates cytosolic channellike structures in which diffusion is accelerated because ions cannot spread in all directions.

Discussion
IP 3 plays a major role in Ca 2+ signaling by mobilizing Ca 2+ from the ER, which is the main intracellular Ca 2+ store. After stimulation, IP 3 must diffuse from the plasma membrane where it is produced across the cytoplasm to trigger Ca 2+ release from the ER. Thus, IP 3 diffusion plays an essential role in Ca 2+ signaling. Following measurements in a medium devoid of IP 3 R, IP 3 diffusion was assumed to be fast, with a diffusion coefficient 10 of 283 ± 53 μm 2 s −1 . This estimation is in line with those of compounds of similar molecular weight such as ATP for example 40 . However, this value cannot account for the observation that in response to a localized release of a nonmetabolizable analog of IP 3 , there is a ~ 10 s delay in the occurrence of Ca 2+ puffs 30 μm away from the site of IP 3 release 11 . In the same manner, such a rapid diffusion cannot account for the drastic influence of the localization of the IP 3 -metabolizing 5-phosphatase enzyme on Ca 2+ signaling observed in CHO cells by De Smedt et al. 41 . The latter authors developed a mutant InsP 3 5-phosphatase in which the C-terminal cysteine cannot be farnesylated, which hinders its binding to the plasma membrane. While the Ca 2+ oscillations detected in the presence of 1 μM ATP were totally lost in 87.5% of intact (farnesylated) InsP 3 5-phosphatase-transfected cells, a loss of Ca 2+ signal occurred in only 1.1% of the mutant InsP 3 5-phosphatase-transfected cells 41 . Such a sensitivity to the location of  A and B) show the 2D geometry considered in the computational simulations of IP 3 diffusion and Ca 2+ puff occurrence in response to the localized photorelease of a non-metabolizable IP 3 analogue. The shape of the cell was redrawn from Dickinson et al. (2016) and that of the ER was largely inspired from De Angelis et al. 39 . Panels (C and D) show the simulated (plain symbols) and experimental (empty symbols) mean latencies (blue) and minimal first puff latencies (black). Flash duration is 500 ms. Simulation procedures are the same as for Fig. 4. The rates of localized IP 3 photorelease, θ, is taken equal to 400.6 μMs −1 and 284.65 μMs −1 , respectively which corresponds to the 250 μMs −1 used in the other simulations with the volumetric adjustments. www.nature.com/scientificreports/ an IP 3 metabolizing enzyme could not be observed if IP 3 was indeed such a fast-diffusing molecule that rapidly becomes homogeneously distributed in the whole cell. Accordingly, a much lower value of the effective diffusion coefficient of IP 3 was predicted by Dickinson et al. 11 on the basis of their observations of puff latencies in SH-SY5Y cells in response to the photorelease of a caged non metabolisable IP 3 analog. Although this indirect evaluation of the rate of IP 3 spreading is submitted to uncertainties due to the stochastic nature of Ca 2+ puffs and to their dependence on Ca 2+ concentration that is not fully controlled by EGTA, it has the unique advantage of corresponding to realistic cellular conditions. Based on a simplified description of Ca 2+ puff dynamics, these authors proposed that observations can be accounted for if the effective diffusion coefficient of IP 3 in these cells is in the 3-10 μm 2 s −1 range. This conclusion was reexamined in the present study, based on a detailed description of Ca 2+ puff dynamics. Our results also point to a lower value of the effective diffusion coefficient of IP 3 in intact cells than in cytosolic extracts of Xenopus oocytes but concluded to a decrease by a factor ~ 3 instead of ~ 30. On the basis of their in vitro measurement, Allbritton et al. 10 predicted the range of action of IP 3 to be ~ 24 μm. This is an estimation of the range on which IP 3 can diffuse before being metabolized into another inositol phosphate. Given that this range is proportional to √ D I , we would predict a range of action of 14 μm. This value is estimated at 4.5 μm with the D I value inferred in Dickinson et al. 11 . Since most mammalian cells have a diameter of ~ 12 μm 42 , the two values lead to significantly different physiological conclusions. With an effective diffusion coefficient of 100 μm 2 s −1 , IP 3 acts as a global messenger in many cases, although this would not be the case in very large cells such as oocytes or some adipocytes.
Besides the different computational framework, our approach differs from that developed by Dickinson et al. (2016) in three ways. First, only one spatial dimension was considered in the latter study, while we performed 2D or 3D simulations. As illustrated for simulations of IP 3 diffusion in the 1D-like astrocytic process ( Fig. 7 and Video S5), diffusion is faster in 1D. The 1D approach of the previous study led to an underestimation of the diffusion coefficient because the value was fitted to reproduce the effective rate of IP 3 diffusion that was observed in 3D. A second important difference relates to the relation between IP 3 concentration and the probability of Ca 2+ puff occurrence. At low IP 3 concentration, the linear relation considered in the previous study predicts a www.nature.com/scientificreports/ larger probability of puff occurrence than the nonlinear function used in the present study (Eq. 2). Again, this higher probability of puff occurrence was compensated by a lower value of the IP 3 diffusion coefficient. Thirdly, because we here considered both the mean and the minimal first puff latencies, we inferred the rate of IP 3 release by the flash (parameter θ in Eqs. 3 and 5) in a more accurate way. We reasoned that the mean first puff latency is affected by an increase in Ca 2+ because the latter is not fully prevented by the EGTA injected in the cell at a final concentration of 5 μM. In line with this hypothesis, agreement for both mean and first puff latencies was only obtained in our simulations when considering that locally the level of Ca 2+ was in average 120 nM above basal level when puff activity was monitored after distributed photorelease of IP 3 . Because this Ca 2+ increase was not taken into account in the previous study of Dickinson et al. 11 , the value of θ was overestimated, leading again to an underestimation of the IP 3 diffusion coefficient.
To confirm that the IP 3 diffusion coefficient is indeed larger than 10 μm 2 s −1 , we simulated the experimental protocol of distributed photorelease of IP 3 used in Dickinson et al. 11 to induce a spatially uniform rise in IP 3 concentration and described in the methodology section here above. Simulations performed with D I = 10 μm 2 s −1 and shown in Figure S6B&D show that IP 3 concentration remains spatially inhomogeneous up to at least 10 s after the flash. This is not in agreement with the observation that first puff latencies are independent from the distance from cell end under this protocol. In contrast, when D I = 100 μm 2 s −1 , IP 3 rapidly equilibrates in the whole cell ( Figure S6C).
Simulations predict that IP 3 spreading is not much affected by the proximity of the plasma membrane. Indeed, times of puffs occurrence are mainly determined by the distance of the clusters from the spot of IP 3 photorelease and not by their positions in the bulk of the cytoplasm (as in Fig. 4) or close to the plasma membrane (as in Fig. 6). This is relevant because clusters of IP 3 receptors are most of the time located very close to the plasma membrane 37 . Similarly, IP 3 spreading is not much affected by the precise shape of the cell. Simulations of puff activity indeed lead to similar results in a 2D ellipse-shaped cell (Fig. 4), in a 3D ellipsoid (Fig. 6) or in a 2D configuration resembling a SH-SY5Y cell (Fig. 8). However, elongated geometries favor fast diffusion and allow for higher local [IP 3 ] than in extended systems. This was shown here by simulations in a situation corresponding to an astrocyte. Accordingly, rates of Ca 2+ waves propagation in astrocytic processes are larger than in most cell types (33-100 μms −1 , see Cornell-Bell and Finkbeiner 43 versus 10-20 μms −1 , see Dupont et al. 23 ). Surprisingly, the presence of the ER, which acts as an obstacle to free diffusion, does not much affect mean effective diffusion. This lack of global effect can be ascribed to the fact that while some paths are slower because of obstruction, others are faster because of the presence of channel-like structures which favor fast diffusion thanks to the absence of dilution effect. It should be noted that we did not consider the tubular-shaped of this network, nor its dynamical evolution into a more fragmented structure 44 .
The value of 100 μm 2 s −1 for the effective diffusion coefficient of IP 3 that emerges from our simulations is in line with the buffering effect exerted by the non-fully bound IP 3 R tetramers, as proposed previously 11,13 . Thus, it is expected to be smaller in cell types with higher expression levels of IP 3 R. In the same line, effective diffusion is much accelerated upon increasing [IP 3 ] because buffers become saturated. Thus, upon cell stimulation by an agonist that generally leads to a surge in [IP 3 ] followed by a decrease, the properties of diffusion are expected to vary. Even more significant in this respect are the IP 3 oscillations that have been observed in several cell types and that arise from the activation by Ca 2+ of IP 3 synthesis by PLC 45 and/or of IP 3 metabolism by a 3-kinase 46,47 . Further studies are required to assess the consequences of the interplay between the temporal changes in IP 3 concentration and its diffusional properties.